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Collocation 
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John W. Dankanich 
AeroDank, Inc. 

Cleveland, Ohio 44135 


As NASA examines potential missions in the post space shuttle era, there has been 
a renewed interest in low-thrust electric propulsion for both crewed and uncrewed mis- 
sions. While much progress has been made in the field of software for the optimization 
of low-thrust trajectories, many of the tools utilize higher-fidelity methods which, while 
excellent, result in extremely high run-times and/or poor convergence when dealing with 
planetocentric spiraling trajectories deep within a gravity well. Conversely, faster tools 
like SEPSPOT provide a reasonable solution but typically fail to account for other forces 
such as third-body gravitation, aerodynamic drag, solar radiation pressure. SEPSPOT is 
further constrained by its solution method, which may require a very good guess to yield 
a converged optimal solution. Here the authors have developed an approach using collo- 
cation intended to provide solution times comparable to those given by SEPSPOT while 
allowing for greater robustness and extensible force models. 


Nomenclature 


t Time, days 

x State variable vector 

u Time- varying control variable vector 

(ft g Global design variable vector 

<ft p Phase design variable vector 

ifto Initial boundary constraint 

iftf Final boundary constraint 

xftp Path constraint 

r Non-dimensional time within a segment polynomial 

Ac Continuity defect 

As Differential defect 

Li Lagrange interpolation matrix for interior nodes 

D, Lagrange derivative matrix for interior nodes 

D c Lagrange derivative matrix for cardinal nodes 

n Number of state variables 

m Number of time-varying control variables 

ncn p Number of cardinal nodes per phase 

nirip Number of interior nodes per phase 

z Independent variable vector in the nonlinear optimization problem 

g Boundary constraint vector for the nonlinear optimization problem 

a Semi-major axis, km 

e Eccentricity 

i Inclination, deg 
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Q 


Right ascension of ascending node, deg 
u> Argument of periapsis, deg 

8 True anomaly, deg 

h Specific angular momentum, ^ L - 

p Semilatus rectum, km 

/ X-component of eccentricity vector in the equinoctial frame 

g Y-component of eccentricity vector in the equinoctial frame 

j X-component of ascending node vector in the equinoctial frame 

k Y-component of ascending node vector in the equinoctial frame 

L True longitude, deg 

r Spacecraft radius from central body, km 

v Spacecraft velocity in the EME2000 frame, 

m Spacecraft mass, kg 

rev Number of orbital revolutions 

T p Orbital period, days 

PQW Perifocal coordinate frame 

NTW Velocity-vector aligned spacecraft coordinate frame 
RSW Local-vertical, local-horizontal spacecraft coordinate frame 

EME 2000 Earth Mean Equator of J2000 inertial coordinate frame 
r c t Central body equatorial radius, km 

p Gravitational parameter of the central body, 

Po Spacecraft propulsion system input power, kW 

77 Spacecraft propulsion system efficiency 

I sp Spacecraft propulsion system specific impulse, sec 

Subscript 

L Lower bound 

U Upper bound 

i Interior node 

c Cardinal node 

0 Initial condition 

f Final condition 


I. Introduction 

F OR several decades, SEPSPOT has been NASA’s primary tool for the optimization of planetocentric 
low-thrust trajectories. 1 SEPSPOT takes advantage of orbital-averaging techniques to greatly increase 
the speed of computation at the expense of fidelity. Generally, orbital averaging techniques provide an 
accurate answer regarding the transfer of one orbit to another, but lack the precision to put a spacecraft at 
a given position at a specific time. Although relatively fast, SEPSPOT suffers from some weaknesses that 
result in difficulty when using it to analyze state-of-the-art low-thrust transfers. SEPSPOT utilizes indirect 
optimization techniques based on Hamilton-Lagrange theory. In addition to providing an initial guess at 
values of the orbital state variables, one must also provide a fairly accurate guess for initial and final values 
of the costates, which can be challenging. More critically, extending the equations of motion to include 
additional force models (atmospheric drag, solar pressure, perturbations by the Moon and Sun, etc.) would 
require the equations of motion of the costate variables to be rederived from their current form. 

To address these issues, the authors have embarked on developing an approach which uses a collocation 
technique to transform an initial- value-problem (IVP) or two-point boundary value problem (TPBVP) into a 
nonlinear programming (NLP) problem, which may be solved utilizing an off-the-shelf optimization software 
package. Since this approach utilizes direct optimization method, it has the advantage that new environ- 
mental and vehicle related force models may be added with relative ease. Collocation and pseudospectral 
optimization techniques were demonstrated by Dickmanns and Well 2 and Hargraves and Paris 3 and have 
been used with success in software such as OTIS4, 4 DIDO, 5 and SOCS. 6 Collocation and pseudospectral 
methods have previously been applied to the low-thrust spiraling problem by Betts.' However, in his work 
he directly transcribed the equinoctial equations of motion with respect to time, which led to a mesh grid 
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of over 16000 points and a run time of several hours. 7 Rather than take this approach, which sees highly 
oscillatory behavior in the time-histories of the states and controls, the authors utilize orbital averaging 
techniques in conjunction with hybrid control formulations which largely remove the oscillatory behavior 
from the problem. This substantially reduces the number and order of the collocating polynomials, which 
reduces the size of the NLP problem and results in run times on the order of a few seconds to minutes. 

Kluever demonstrated an orbital averaging problem which utilizes the Gauss form of the Lagrange plan- 
etary equations in conjunction with a hybrid control technique. 8 Using only a handful of mesh points for the 
controls, he demonstrated the ability to optimize a LEO to GEO transfer in good agreement with SEPSPOT. 
For this approach the authors elected to use collocation of the states and controls since the method tends 
to show an ability to find optimal solutions given a relatively poor initial guess. The authors use equations 
of motion based on a modified set of equinoctial elements, which don’t suffer from as many singularities as 
the Lagrange planetary equations. 

Like Kluever, the authors utilize a hybrid set of control variables. These control variables are the costates 
of the Hamilton-Lagrange formulation. They show a relatively smooth time-history throughout an orbital 
transfer, which is beneficial for collocation, but internally are transformed into a set of pitch and yaw angles 
which may oscillate greatly over the course of any single orbit. Furthermore, due to the extensible nature 
of the collocated form of this problem, alternative guidance strategies can be easily implemented so long as 
they can be parameterized such that their control variables do not exhibit ’’excessively” oscillatory behavior 
throughout the orbital transfer. As such, this method allows users to determine both the optimal (minimum 
time, minimum propellant, etc.) transfer guidance history, but also how it compares to a custom guidance 
strategy. The use of collocation to solve the optimal control problem should, in general, result in a greater 
radius of convergence such that the initial guesses to the costates need not be as accurate as those supplied 
to SEPSPOT. 

This approach does sacrifice fidelity for speed, and thus is intended for the preliminary analysis of low- 
thrust trajectories and guidance algorithm development. For example, orbital averaging is incapable of 
analyzing a lunar insertion maneuver as demonstrated by Betts. However, the run-time is so small that 
using a collocated averaging technique to ’’dig” out of the deepest part of the gravity well would greatly 
reduce analysis time. Once at a sufficiently high intermediate orbit, the oscillatory behavior of the state 
and control time histories is low enough to permit collocation of non-averaged equations of motion from the 
intermediate orbit to the final target. 


II. Formulation 


A. The Collocation Problem 

Collocation is used to simultaneously simulate and optimize the trajectory of a dynamical system. As 
a fallback, the implementation also supports a more traditional explicitly integrated shooting method for 
trajectory optimization, though experience has shown that implicit integration is not only faster, but often 
capable of finding more optimal solutions. 

In general terms, the optimization problem is of the following form: 


Minimize /o&j(x(t 0 | f ), t 0 \ f , u(t 0 | f ), 0 g , (j) v ) (1) 

s.t. x L < x < Xu (2) 

t L <t<tu (3) 

u L < u < uu (4) 

(/>g.L ^ 0g,U (5) 

^p.L 0p f 7 0p,U (6) 

V’o.l < /v> 0 (x(t 0 ),to,u(t 0 ),^g,^p) < ipo,u (7) 

V’f.L < ibo (x(t f ) , tf , u(t f ) ,(f> g ,<f> p ) < V’f.u (8) 

^p.L ^ /pat/i(x(t) , t, u(t) , (^>g , (/>p ) ^ V’p.U (9) 


Where x and u are the time-varying states and controls, respectively. The vectors </> g and 4> p contain 
static optimization variables (design variables) for the entire problem, and a specific phase, respectively. 
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Variable t is the time at the point of the objective function evaluation; either the beginning or end of a 
phase. All states, controls, design variables, and time may have simple bounds. Furthermore, we may 
impose boundary constraints at the start or end of a phase (j/jq and if)/), or path constraints to be assessed 
throughout a phase (pl> p ). 



Time (s) 


Figure 1. The state time-history of a single state in a phase consisting of three 3rd-order segments, including 
the differential defects ( A ,s v ; ) and continuity defects ( Ar v ) 

The collocation problem consists of a series of phases in which the forces acting upon the spacecraft 
are consistent such that the states are C 1 continuous within a phase. Each phase is defined on an interval 
[to, to + t p \ , and the equations of motion dictate some set of states x and controls u. Each phase is further 
subdivided into a series of segments along which the states and controls are represented by a polynomial of 
at least order 3, and thus C 2 continuous. For each state and control along each segment, this polynomial is 
defined at a series of nodes called the cardinal nodes. Internally, each segment of order o is defined on the 
interval [—1, 1] (so called r-space) and has cardinal nodes at the (o + l) Legendre-Gauss-Lobatto points. For 
example, a 3rd-order polynomial segment has four cardinal nodes. The collocation engine uses a nonlinear 
optimization routine to vary the cardinal values of the states and controls, the phase start time and duration, 
and design variable values such that state time histories accurately reflect the equations of motion and the 
objective function (1) is satisfied. 

The accuracy of the collocated polynomials in representing the dynamics of the system is measured by 
comparing the slope of the polynomial for a given state to the derivative of that state as given by the 
equations of motion. These differences are called differential defects. A well-posed initial-value-problem has 
a unique solution and requires that the number of independent variables is equal to the number of equality 
constraints. Figure 1 illustrates the setup of a collocation problem for one state variable in a phase of three 
3rd-order polynomial segments. With just one state this problem has 14 variables: the initial time and phase 
duration (2), and the state values at four cardinal nodes in each of three segments (12). We require that 
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there be no discontinuity in the value of a state variable at the segment boundaries. 


Acij — x cJ x c i — 0 (10) 

Furthermore, the initial time, duration, and initial state value in this phase are fixed via bounds, giving 
a total of five equality constraints. The remaining nine equality constraints are achieved by requiring the 
differential defects at three points in each segment to be zero. The points in each polynomial where the 
defects are assessed are referred to as the interior nodes. Within each segment we may construct a Lagrange 
interpolation matrix Li and a Lagrange derivative matrix D, such that the values and derivatives of the 
states and controls at the interior nodes may be obtained by simple matrix-vector products: 


Xi 

— LiX c 

(11) 

dx j 

d,r 

(12) 

dt 

■ 1 •+■* 
! 1^3 
u 
X 

Q 

II 



(13) 


Finally, the differential defects of the states at the interior node of the phase are: 


AX' 

feom (xi i ti , Uj , 0g , (j) p ) ^ — 0 (14) 

Some preliminary tests suggest that a good choice of interior nodes is to determine the LGL points for 
an (o + 1) polynomial and remove the endpoints. Reusing the first (or last) o cardinal nodes as the interior 
nodes is also possible, but generally exhibits poorer convergence. On the other hand, if the convergence 
issues can be fixed, reusing the cardinal nodes as interior nodes would reduce equation of motion evaluations 
by nearly a factor of two, and no matrix multiplications (11) would be required to determine the state and 
control values at the interior nodes. 

Experience has also shown that convergence is aided by requiring both value and rate continuity to 
be imposed on the control values at the segment boundaries. Control rates at the segment boundaries 
are determined with a Lagrange derivative matrix constructed such that it returns the derivatives of the 
polynomial at the cardinal nodes: 


du c 

dt. 


dr 

D c u c —— 
dt 


(15) 


B. The NLP Interface 

The optimization problem above is transcribed by the collocation engine into the following form: 


Minimize ,fobj( z) (16) 

s-t. g L < g(z) < gu (17) 

z L < z < zu (18) 

This transformed optimization problem may be solved by a variety of off-the-shelf nonlinear optimization 
software packages. Currently the authors are using the IPOPT 9 which is able to capitalize on the sparsity 
of the Jacobian matrix that is characteristic of the collocaiton problem. 

The times, states, controls, and design parameters are scaled and packed into the independent variable 
array for the NLP (z). Currently a basic scaling scheme is used whereby each independent variable is scaled 
by the inverse of its expected magnitude such that the corresponding component of z is roughly on the order 
of one. The independent parameter array is packed in the following order: 

1. The scaled global design parameter values comprise the first elements of z 

2. The scaled design parameter values of the first phase comprise the next elements of z 
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3. The scaled initial time and duration of the first phase are the next two elements of z 

4. The cardinal values of the n states and m controls of the first phase are packed as the next ( n+m)-ncn p 
elements, where ncn p is the number of cardinal nodes in the phase. These values are packed such that 
the states and columns at any given cardinal node are contiguous. Conceptually, this is the same as 
constructing the following ncn p by n + m matrix X C U C and unraveling its values in row-major order. 

* 1,0 * n — l ,0 ^ 0,0 ••• Mra - 1,0 

*1,1 ••• *71-1,1 Wo l ••• W m _l 1 

: . . . . (19) 

\j^0,ncn p — l %l,ncn p — 1 * * * %n— l,ncn p — 1 ^0,ncn p — 1 ' ’ * ^m— l,ncn p — 1_ 

5. Items 2 through 4 are repeated for each phase in the problem. 


z=[l<fi s ] T ttpf to t p [X c u c ] 0 ... [X c U c ] ncnp _J (20) 

Similarly, the constraints of the collocation problem are packed into the NLP constraint array (g) array 
in the following order: 

1. The phase linkage constraints comprise the first elements of g 

2. The nin p ■ n state defects for the first phase are stored node- by- node as the next elements of g. 

3. The state and control continuity conditions at each of the numseg — 1 segment boundaries in the first 
phase make up the next ( numseg — 1) • (n + m) elements of g. 

4. The control continuity rate conditions at each of the numseg — 1 segment boundaries in the first phase 
make up the next ( numseg — 1) • m elements of g. 

5. The initial and final boundary constraints of the first phase make up the next elements of g 

6. Finally, the path constraints are evaluated at each cardinal node and comprise the final ncn p ■ npc 
elements of g, where ncp is the number of path constraints in the phase. 

7. Items 2 through 6 are repeated for each phase in the problem. 


S [^®]():7U7ip — 1 \^'^'\o:num_seg—l \.^^\o:numseg— 1 [M T [V’f ] T [V’pio ••• hML P -i] ( 21 ) 

By ordering the z and g arrays in such a way the sparsity pattern of the Jacobian matrix is largely 
block-diagonal. The use of a nonlinear programming routine which is able to capitalize on the sparsity of 
the Jacobian matrix is a critical feature which enables the solution to the collocation problem to be obtained 
rapidly. 

C. Equations of Motion 

In his work, Kluever showed that, by applying orbital averaging to the Gauss form of the Lagrange planetary 
equations, 10 he was able to get very good agreement with SEPSPOT results by parameterizing the controls 
at some set of control nodes and explicitly integrating the solution. 8 However, the Lagrange planetary 
equations are problematic in the context of a collocation formulation since they exhibit singularities at zero 
eccentricity and zero inclination, a point which is of interest for the purpose of examining trajectories to or 
from geostationary orbit. Instead, the authors utilized a set of modified equinoctial elements as given by 
Walker. 11 Betts also utilized this set of equations of motion in his low-tlrrust orbit transfer example, though 
he did not employ an averaging technique. 12 
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X = 


0 

fsrn (L) 
^cos (L) 
0 
0 
0 

0 

0 


2 £ 


£ [(m+l)cos(L)+/] 
p w 

£ [Q+l)sin(.L)+g] 
fl w 

0 
0 
0 


0 

£ g(jsin(L) — fccos(L)) 
/X w 

£ /(jsin(L)-fccos(L)) 
/I w 

(L) 




/i 2 to k 

p jsin(L) — kcos(L) 
p w 

o 

0 




—m 

rev 


(22) 


where the state vector consists of the modified equinoctial elements and spacecraft mass 


1 T 


p f g j k L m rev 


(23) 


the accelerations 


(24) 


are the radial, local-horizontal, and orbit normal perturbing accelerations where 


s 2 = l+j 2 + k 2 (25) 

w =- = 1 + /cos (L) + gsin (L) (26) 

Note we refer to the equinoctial element h as j to avoid confusion with the specific angular momentum 
of the orbit. The equations of motion (22) account for central-body gravitation in the absence of perturbing 
accelerations. Adding different perturbing accelerations is as simple as converting them to the RSW frame 
and including them in (24). For the purposes of this paper, the only non-two-body accelerations are due to 
thrust and J 2 perturbations. 

The mass flow rate (m) is based on the propulsion system model (29). 

The number of orbital revolutions performed by the spacecraft can be approximated by treating it as an 
integrated state variable with a derivative function equal to the frequency of the spacecraft orbit. 


rev 


1 

% 


27t a 5 


(27) 


D. Solar Electric Propulsion Model 

The solar electric propulsion system generates power based on its distance to the Sun and the angle between 
the array plane and the Sun vector. In the analysis below the Sun is assumed to be 1 A.U. from the spacecraft 
at all times and the arrays are always perfectly pointed at the Sun. If the nominal array power at 1 A.U. is 
given by Pq , then the spacecraft acceleration due to thrust is 8 


CLt 


2rjPo_ 

mgl sp 


and the rate of propellant expenditure due to engine firing is derived from the definition of thrust: 


(28) 


Ft 


ring I> 


sp 


2??Po 

9 2 I 2 p 


(29) 
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E. Gravitational Perturbations due to Oblateness of the Central Body 

Kechichian quantified Ji perturbations in the RSW frame as: 13 


3 ^ 2 ^ ^ _ 12 (j'sin(L)-fccos(L)) 2 ^ 
12p,J 2 r 2 ob (jsin(L) — kcos(L)) (jcos(L)-i-ksin(L)) 

r 4 g 4 

(jsin(L)-fccos(L))(l-fc 2 -j 2 ) 

r- 4 »4 


(30) 


F. Orbital Averaging 

The equations of motion (22) are subject to oscillatory behavior which can greatly increase the number of 
collocation segments required for a converged solution. To enable better performance from the collocation 
routine, the equations of motion are averaged. The averaged equations of motion are obtained by converting 
x to integrating for an entire orbital revolution, and dividing by the orbital period to obtain the average 
rates of change in the modified equinoctial elements for a given orbital state. 8 The anomaly term (true 
longitude for the modified equinoctial elements) is not included among the averaged orbital elements. 


x = 


i_ r 

T p dt dd 


(31) 


The transformation from time to true anomaly ( 6 ) is approximated using the perturbation-independent term 
of the Gauss form of the Lagrange planetary equation for true anomaly. 10 


d9 _ h_ 
dt r 2 


(32) 


G. Eclipse Arcs 

While accelerations due to solar electric propulsion only need to compute the integral in (31) through the 
illuminated orbital arcs, forces due to gravitational harmonics, solar radiation pressure, drag, third-body 
effects, and trapped radiation impingement require that the entire orbit be integrated. To determine those 
points at which the spacecraft enters and exits the shadow of the central body, the authors employed Escobal’s 
quartic shadow function: 1415 


S = opcos (9) + ccicos (9) 3 + c^cos (9) 2 + oqcos (9) + oq 


(33) 


where 


Q. 1 


e 4 

-2(/?f - /3i 2 )e 2 

(/3 2 2 + /?r 2 ) 2 

OL 2 


4e 3 

-4e(/3| - /3?) 

0 

Oi 3 

= 

6e 2 

— 2((/?2 ~ Pi) + e 2 * (1 — /3f)) m 

-Pi){i-PD-mp 2 ) 

«4 


4e 

-4e(l - /? 2 ) 

0 

. 0 : 5 . 


. 1 

-2(1 -ft) 

(1-/3I) 2 


a = 


r c b 


(34) 


(35) 


and fi is the unitized vector from the central (shadowing) body to the Sun in the perifocal frame (PQW): 15 

G = PiP + /3 2 q + /33W (36) 
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The solution of the quartic shadow funtion (33) will yield false positives, which must be filtered out. 
Furthermore, since we seek 0 and the solution is in the form cos (9), an even function, we must account for 
all possible combinations. This yields as many as eight potential solutions for the roots, expressed in terms 
of sin (9) and cos (9) : 

COs(0i) COS (02) COS (0 3 ) COS ($4) COS (01 ) COS (02) cos (03) cos (04 ) (37) 

sin(0!) sin(0 2 ) sin(0 3 ) sin(0 4 ) — sin^) — sin(0 2 ) — sin(0 3 ) — sin(0 4 ) (38) 

The true roots in cos (0) are obtained by passing (37) through the shadow function as given by Vallado 
(39). 15 Here the shadow function has been normalized by the central body radius to help with numerical 
errors: 

5(0) = (1 + ecos (0)) 2 + (ft cos (0) + ft sin (0)) 2 - (39) 

Both Escobal’s quartic shadow function (33) and the form given in Vallado (39) will also have valid roots 
on the illuminated side of the central body. These are filtered by assuring that we only take those roots for 
values of true anomaly at which central body is between the spacecraft and the Sun. Using equation 5-3 
from Vallado, 15 the shadow entry and exits must be such that: 


f3 icos (0) + ftsin (0) < 0 (40) 

Now the roots of the shadow function in true anomaly are known, and the derivative of the shadow 
function at the roots indicates whether the root is a shadow entry (^f > 0) or shadow exit (^ < 0). 

^ = 2 (j^j (ft cos (0) + ftsin (0))(— ftsin (0) + ftcos (0)) - 2esin (0) (1 + ecos (0)) (41) 

We can introduce breakpoints in the integral from (31) and evaluate the integral along each illuminated 
or shadowed segment, applying thrust due to a solar electric propulsion system accordingly. For each interval 
in a given orbit, the integrand in (31) is evaluated using a ninth-order Legendre-Gauss quadrature. 

H. Control Parameterization 

The forces due to thrust are determined by the guidance scheme of the spacecraft. In his non-averaged 
scheme, Betts used a unitized thrust vector in the RSW frame as the control. 12 This is not possible using 
an orbital-averaging approach, since a single control vector is needed that will provide the appropriate 
oscillatory behavior of the thrust direction within a single averaged orbit. Kluever showed that the costates 
of the classical orbital elements can be used as the time-varying guidance parameters for an averaging 
approach. 8 Kluever’s approach (42) assumes we only have costates governing the ’’slow” orbital elements, 
and provides a thrust vector in the velocity vector-aligned frame (NTW): 


0 2 — 

rsin(0) 2(e+cos(0)) 

a-v v 

arMcoefl = o 0 

0 0 

2e +( r £ 2 ^g)) 2sin(g) 

e-v e-v 


o 

o 

cosOy+S) 

' h 

sin(aj+fl) 
hsin(inc ) 

— rsin(u;+0)cos(mc) 
hsin(inc ) 



(42) 
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Since our equations of motion are formulated using modified equinoctial elements, the classical elements 
must first be computed (see reference 12). The guidance parameters are a given as a vector of the classical 
orbital element costates, excluding true anomaly. 


^coe — 




A a X e X i Xq A u 


and the thrust unit vector in the NTW frame is: 


l NTW 


M-coe ^coe 

l|Mj oe A coe || 


(43) 


(44) 


In his work, Kluever divides A 0 by the semi-major axis and parameterizes the other costates as a function 
of the semi-major axis. 8 In this approach A a is normalized by the semi-major axis before being used in (44), 
but all costates are functions of time. Notably, since the guidance parameters here are not functions of 
semi-major axis, the posibility exists to use these controls for a single-phase ’’round-trip” trajectory or other 
trajectories which involve both increasing and decreasing the semi-major axis of the orbit. 

We may derive a similar guidance scheme where the equinoctial equations of motion and costates are 
used in place of the classical orbital element equations and costates. 


®TM mee ii — 


£sin (L) 


^cos (L) 


■2p p 

w y p 
p [(-i«+l)cos(L)+/] 


p [(w+l)sin(L)+g] 


p g(jsin(£)-fccos(L)) 
fl w 

p /(jsin(L)-fccos(L)) 


^cos(L) 




(45) 


The guidance parameters are a given as a vector of the modified equinoctial element costates, excluding true 
longitude. 


A 


mee 


A p 


•V ^ g 



and the thrust unit vector in the RSW frame is: 


(46) 


— 



(47) 


In the modified equinoctial element control formulation, X p is divided by the senri-latus rectum before 
being passed to (47). 

The modified equinoctial element-based guidance scheme is similar to that used in SEPSPOT 1 and is 
less prone to singularities, except when an orbit is nearly retrograde. However, if one only cares to control 
semi-major axis, eccentricity, and inclination of the spacecraft, this requires only three time-varying control 
variables (A a , A e , A;) in the case of the classical element formulation, but five time-varying control variables 
(Xp, A/, X g , Xj, Afc) in the case of of the modified equinoctial element formulation. The additional control 
variables increase the size of the collocation problem, and may adversely affect runtime, though this effect 
has not yet been quantified. 


III. Results 

The collocation formulation of a planetocentric low-thrust trajectory optimization problem shown here 
is demonstrated with comparisons to the LEO to GEO and GTO to GEO cases in Kluever’s demonstration 
of his solution method. 8 
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A. Minimum Time LEO to GEO 


The minimum time LEO to GEO spiral transfer uses the same assumptions as those given by Kluever. 8 The 
spacecraft initially has a mass of 1200 kg and the propulsion system has a nominal power level of 10 kW, 
a specific impulse of 3300 sec, and a propulsive efficiency of 65%. The initial time of the propagation is 
January 1, 2000 and the initial state vector is given in Table 1. 


Variable 

Initial Value 

P 

6927 km 

f 

1.0E-6 

9 

0 

j 

0.2539676 

k 

0 

m 

1200 kg 

rev 

0 


Table 1. LEO to GEO initial variable values (fixed) 


The solution shown here uses the classical orbital element costates as guidance parameters. The semi- 
major axis costate (A a a) is fixed at -1, giving a tangential thrust vector when A e = 0 and A = 0. Costates 
for right ascension of ascending node (An) and argument of periapsis (A w ) are fixed at 0. 

The costates corresponding to eccentricity and inclination are time-varying optimal controls. Initially 
the value of these controls is a linear fit of the guessed initial and final values. 


Parameter 

Initial Value (guess) 

Final Value (guess) 

Lower Bound 

Upper Bound 

Scale Factor 

A e 

0 

0.5 

-1 

1 

100 

Ai 

0 

3 

-10 

10 

100 


Table 2. LEO to GEO time varying guidance parameters 


Initially, the trajectory is simulated starting at the initial state for the guessed elapsed time using linear 
fits to the time-varying controls. This explicit simulation provides values for the state variables at the 
cardinal nodes in the phase. At this point, the solution is physically accurate (the defects are approximately 
zero) but the constraints are not yet satisfied. 


Variable 

Lower Bound 

Upper Bound 

Scale Factor 

a 

42164 (km) 

42164 (km) 

1 

42164 

e 

0.0001 

0.001 

100 

i 

0 (deg) 

0.01 (deg) 

100 


Table 3. LEO to GEO final boundary constraints 


In addition to the boundary constraints above, a path constraint is used to ensure the value of periapsis 
altitude is at least 300 km. Without this path constraint, the optimizer sometimes attempts to push periapsis 
of the orbit below the surface of the Earth. 

With the initial trajectory of the vehicle reasonably defined from the explicit simulation and the boundary 
and path constraints in place, the solution is solved using IPOPT. 9 Reasonable convergence was achieved 
with the phase broken into 15 equal 3rd-order polynomial segments. The solution converged in 91 iterations. 
The resulting minimum transfer time is given in Table 4 including comparisons to results obtained by Kluever 
and SEPSPOT. 

The time histories of semi-major axis, eccentricity, and inclination are shown in figure 2. They show good 
agreement with the results achieved by Kluever. 8 Figure 3 shows the time history for the eccentricity and 
inclination costates used as optimal control variables. Both exhibit some ’’wagging” at the beginning and end 
of the phase. This indicates that the scaling in the collocation problem needs improvement, but the very good 
agreement with other results suggests that the solution is insensitive to the values of the control variables at 
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those times. Better scaling may reduce or eliminate the ’’wagging” behavior and improve convergence but 
not have a significant impact on the solution. 




Figure 2. Semi-major axis, eccentricity, and inclination histories for the minimum-time LEO to GEO transfer. 
Markers indicate cardinal nodes of the collocation problem. 



Figure 3. Eccentricity and inclination costates used as the time-varying optimal controls for the minimum-time 
LEO to GEO transfer. 


Solution Technique 

Transfer Time (days) 

Collocation 

198.6 

Kluever Control Parameterization 8 

199.0 

SEPSPOT Result 8 

198.8 


Table 4. LEO to GEO results 
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Figure 4 shows the evolution of the orbit from LEO to GEO for the minimum time solution, with orbits 
plotted at the collocation segment boundaries, roughly once every 13.2 days. 






Figure 4. Evolution from the initial orbit (blue) to the final orbit (red) for the minimum time LEO to GEO 
transfer. 


B. Minimum Time GTO to GEO 

As in the previous case, the minimum time GTO to GEO spiral transfer uses the same assumptions as those 
given by Kluever. 8 The spacecraft initially has a mass of 1200 kg and the propulsion system has a nominal 
power level of 5 kW, a specific impulse of 1800 sec, and a propulsion system efficiency of 55%. The initial 
time of the propagation is March 22, 2000 and the initial state vector is fixed to the values given in Table 5. 

The costates corresponding to eccentricity and inclination are time-varying optimal controls with their 
values linearly fit to guessed initial and final values (Table 6), and the initial state values at the cardinal 
nodes are obtained through an explicit simulation. The final boundary constraints are the same as those 
given in Table 3 and the path constraint on periapsis radius has the minimum value constrained to 185 km. 

In this case the phase was broken into 10 equal 3rd-order polynomial segments. Despite the reduction 
in the number of variables in the problem due to the reduced number of cardinal nodes, the GTO to GEO 
transfer required 299 iterations to achieve convergence. The resulting minimum transfer time is given in 
Table 7 with comparisons to results obtained by Kluever and results generated by SEPSPOT. 
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Variable 

Initial Value 

P 

11359.07 km 

f 

0.7306 

a 

0 

j 

0.2539676 

k 

0 

m 

1200 kg 

rev 

0 


Table 5. GTO to GEO initial time and state values (fixed) 


Parameter 

Initial Value (guess) 

Final Value (guess) 

Lower Bound 

Upper Bound 

Scale Factor 

A e 

0.5 

6 

0 

15 

100 

Aj 

1 

10 

-15 

15 

100 


Table 6. GTO to GEO time-varying guidance parameters 


Solution Technique 

Transfer Time (days) 

Collocation 

118.29 

Kluever Control Parameterization 8 

118.36 

SEPSPOT 

118.29 


Table 7. GTO to GEO results 


Figure 5 shows the time histories of semi-major axis, eccentricity, and inclination for the minimum time 
GTO to GEO solution, which again compare favorably with the results achieved by Kluever. 8 




Figure 5. Semi-major axis, eccentricity, and inclination history for the minimum-time GTO to GEO transfer. 


Figure 6 shows the time history for the eccentricity and inclination costates used as optimal control 
variables. Again, ’’wagging” behavior is present, especially in the inclination costate. The insensitivity of 
the result to the terminal value of the costates is a likely source of some of the convergence issues experienced 
in this case. One possible reason for this is scaling. However, the fact that the behavior seems to be exhibited 
where the orbit is nearly circular or near zero inclination suggest the singularities in the classical orbit element 
control formulation may be responsible. 
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Figure 6. Eccentricity and inclination costates used as the time- varying optimal controls for the minimum-time 
GTO to GEO transfer. 


Figures 7 and 8 show the evolution of the orbit from GTO to GEO for the minimum time solution. Orbits 
are plotted at the collocation segment boundaries, or approximately once every 12 days. 




EME2000 x (km) EME2000 y (km) 

Figure 7. Evolution from the initial orbit (blue) to the final orbit (red) for the minimum time GTO to GEO 
transfer. 
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Figure 8. Evolution from the initial orbit (blue) to the final orbit (red) for the minimum time GTO to GEO 
transfer. 


IV. Conclusions and Future Work 

The results shown here indicate that a collocation technique, combined with orbital averaging of the 
modified equinoctial elements, is capable of generating optimal solutions for planetocentric low-thrust orbit 
transfer problems relatively quickly. While the results given here were set up for comparisons with Kluever’s 
results, 8 the implementation can be extended to include other forces on the spacecraft such as atmospheric 
drag and perturbation by the Moon and Sun. 

Despite the fact that this approach uses costates of the classical or modified equinoctial orbital elements, 
the use of collocation instead of a shooting method generally allows for greater convergence. Rough guesses 
at the guidance parameters can yield an optimal solution. 

The issues of high nonlinearity in the terminal control conditions clearly need to be addressed. Better 
scaling of the NLP will likely help, but experience with OTIS4 has shown that other measures may be 
necessary. In the use of OTIS4, a common approach to such issues is to add control rate constraints, which 
often forces the control time history to be more smooth. A grid refinement algorithm should also be put 
into place to ensure that the converged solution adequately matches an explicitly integrated trajectory. 

The ease of which this method can be extended to include other control formulations lends itself to 
experimenting with more control laws in the future. Solving these problems using the modified equinoctial 
element control formulations may also prove to be more robust due to the lack of singularities in those 
equations. Other guidance laws which do not require high-frequency time- varying controls, such as Q-Law, 16 
will also be explored. 

Modeling the spacecraft subsystem in greater fidelity will also be a priority. Degradation of the solar 
arrays due to the impingement of trapped particles can be implemented as it was in SEPSPOT. Reference 1 
includes a great amount of detail as to that implementation. The authors intend to implement more detailed 
radiation modeling in the future using more recent and higher- fidelity models. 

So far this method has been demonstrated for a single phase trajectory. Using a multiple phase approach 
in which later phases employ non-averaged equations of motion may allow for fast, robust solutions to 
planetocentric spiral trajectories which target lunar orbit or lagrange points in a terminal phase. 
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